library(prettydoc)
library(tidyverse)
library(gsw)
library(readr)
hydrostation_bottle <- read_delim("hydrostation_bottle.txt",
delim = "\t", escape_double = FALSE,
col_names = FALSE, trim_ws = TRUE, skip = 31)
hydrostation_bottle_names <- read_csv("hydrostation_bottle.txt",
skip = 30)
colnames(hydrostation_bottle) = colnames(hydrostation_bottle_names)
#view(hydrostation_bottle)
# lets first plot the data
hydrostation_bottle %>%
filter(`Sig-th`!= -999) %>% # filter out -999 no data flag
ggplot()+geom_point(aes(x=decy, y = `Sig-th`)) # hard to interpret even w no -999s
hydrostation_bottle %>%
filter(`Sig-th`!= -999 & Depth <20) %>% # filter out -999 no data flag and by upper 20m
ggplot()+geom_line(aes(x=decy, y = `Sig-th`)) # line shows seasonality
# clear seasonal signal for sigma-theta, lets see how this compares to temp
hydrostation_bottle %>%
filter(`Sig-th`!= -999 & Depth <20) %>% # filter out -999 no data flag and by upper 20m
ggplot()+geom_point(aes(x=Temp, y = `Sig-th`))
?gsw #launches documentation for gibbs seawater toolbox
## starting httpd help server ... done
? gsw_sigma0 # lets check this function
# it says we need absolute salinity and conservative temperature
#first we need absolute salinity
?gsw_SA_from_SP
#practical salinity
#sea pressure(dbar)
#longitude
#latitude
#plot our pressure data - its missing before 1980s
hydrostation_bottle %>%
ggplot()+geom_point(aes(x=decy,y=Pres))
#we have depth for the time series
hydrostation_bottle %>%
ggplot()+
geom_point(aes(x=decy, y=Depth))
#adds a pressure column from depth and latN columns from/to hydrostation bottle
hydrostation_bottle =
hydrostation_bottle %>%
mutate(Pres_gsw = gsw_p_from_z(Depth*-1,latN,))
hydrostation_bottle %>%
ggplot()+geom_point(aes(x=Pres,y=Pres_gsw))
## Warning: Removed 3 rows containing missing values (`geom_point()`).
# we see strong 1:1 agreement between measured pressure and calculated pressure
hydrostation_bottle %>%
ggplot()+geom_point(aes(x=decy,y=latN))
hydrostation_bottle =
hydrostation_bottle %>%
mutate(Pres_gsw = gsw_p_from_z(Depth*-1,latN,)) %>%
mutate(S_abs_gsw=gsw_SA_from_SP(Sal1,Pres_gsw,360-lonW,latN))
#plot it
hydrostation_bottle %>%
filter(Sal1!=-999) %>%
ggplot()+
geom_point(aes(x=Sal1,y=S_abs_gsw))
# now we need to calculate conservative temperature
# we need absolute salinity, insitu temp (ITS-90), and sea pressure
hydroS =
hydrostation_bottle %>%
filter(Sal1!=-999) %>%
mutate(Pres_gsw = gsw_p_from_z(Depth*-1,latN,)) %>%
mutate(S_abs_gsw=gsw_SA_from_SP(Sal1,Pres_gsw,360-lonW,latN)) %>%
mutate(T_cons_gsw = gsw_CT_from_t(S_abs_gsw,Temp,Pres_gsw))
hydroS %>%
filter(Temp!= -999) %>%
ggplot()+
geom_point(aes(x=Temp, y =T_cons_gsw))
#add line to calculate conservative temperature
hydroS =
hydrostation_bottle %>%
filter(Sal1!=-999) %>%
mutate(Pres_gsw = gsw_p_from_z(Depth*-1,latN,)) %>%
mutate(S_abs_gsw=gsw_SA_from_SP(Sal1,Pres_gsw,360-lonW,latN)) %>%
mutate(T_cons_gsw = gsw_CT_from_t(S_abs_gsw,Temp,Pres_gsw)) %>%
mutate(Sig_th_gsw = gsw_sigma0(S_abs_gsw,T_cons_gsw))
hydroS %>%
filter(`Sig-th`!= -999) %>%
ggplot()+
geom_point(aes(x=`Sig-th`, y=Sig_th_gsw))
hydroS %>%
filter(Sig_th_gsw<0)
## # A tibble: 544 × 19
## Id yyyymmd decy time latN lonW Depth Pres Temp CTD_S Sal1 Sig-t…¹
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 6.00e9 1.96e7 1955. 1703 32.2 64.4 495 -999 -999 -999 36.3 -999
## 2 6.00e9 1.96e7 1955. 1601 32.1 64.3 144 -999 -999 -999 36.6 -999
## 3 6.00e9 1.96e7 1956. 1704 32.1 64.3 873 -999 -999 -999 35.4 -999
## 4 6.00e9 1.96e7 1956. 1704 32.1 64.3 50 -999 -999 -999 36.2 -999
## 5 6.00e9 1.96e7 1956. 1704 32.1 64.3 244 -999 -999 -999 36.5 -999
## 6 6.00e9 1.96e7 1956. 1704 32.1 64.3 685 -999 -999 -999 36.0 -999
## 7 6.00e9 1.96e7 1956. 1705 32.2 64.3 50 -999 -999 -999 36.4 -999
## 8 6.00e9 1.96e7 1957. 1703 32.1 64.3 2024 -999 -999 -999 35.0 -999
## 9 6.00e9 1.96e7 1958. 1504 32.0 64.4 961 -999 -999 -999 35.2 -999
## 10 6.01e9 1.96e7 1958. 1705 32.2 64.2 553 -999 -999 -999 36.1 -999
## # … with 534 more rows, 7 more variables: `O2(1)` <dbl>, OxFix <dbl>,
## # Anom1 <dbl>, Pres_gsw <dbl>, S_abs_gsw <dbl>, T_cons_gsw <dbl>,
## # Sig_th_gsw <dbl>, and abbreviated variable name ¹​`Sig-th`
#view()
# how to replace data in a line without making this many dataframes (i.e. select and/or filters)
hydroS_correctedS_a =
hydroS %>%
filter(Sig_th_gsw<0) %>%
mutate(S_abs_gsw=gsw_SA_from_SP(CTD_S,Pres_gsw,360-lonW,latN)) %>%
mutate(T_cons_gsw = gsw_CT_from_t(S_abs_gsw,Temp,Pres_gsw)) %>%
mutate(Sig_th_gsw = gsw_sigma0(S_abs_gsw,T_cons_gsw))
hydroS_correctedS_b =
hydroS %>%
filter(Sig_th_gsw>0)
hydroS_corrected = rbind(hydroS_correctedS_a, hydroS_correctedS_b)
hydroS_corrected %>%
filter(`Sig-th`!= -999) %>%
ggplot()+
geom_point(aes(x=`Sig-th`, y=Sig_th_gsw))
## Next day
hydroS_corrected %>%
ggplot()+
geom_point(aes(x=Sig_th_gsw, y = Depth))+
scale_y_reverse()+
scale_x_continuous(position ="top")+
xlab(expression(paste(sigma[theta]," (kg m"^"-3",")")))+
ylab("Depth(m)")
## Warning: Removed 543 rows containing missing values (`geom_point()`).
hydroS_shallow = hydroS_corrected %>%
filter(Depth<30)
sig_theta_tm = lm(Sig_th_gsw~decy,data = hydroS_shallow) # y (Sig) is a function of x(decy), lm(y~x,data = data )
summary(sig_theta_tm)
##
## Call:
## lm(formula = Sig_th_gsw ~ decy, data = hydroS_shallow)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5789 -0.8683 0.1070 0.8819 1.5805
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.4047230 1.6704055 19.998 < 2e-16 ***
## decy -0.0042378 0.0008387 -5.053 4.59e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9401 on 3410 degrees of freedom
## (56 observations deleted due to missingness)
## Multiple R-squared: 0.007431, Adjusted R-squared: 0.00714
## F-statistic: 25.53 on 1 and 3410 DF, p-value: 4.587e-07
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
plot = hydroS_shallow %>%
ggplot()+
geom_point(aes(x=decy, y=Sig_th_gsw))+
geom_line(aes(x=decy, y=Sig_th_gsw))+
geom_smooth(aes(x=decy, y=Sig_th_gsw),method = "lm")+
theme_classic()
ggplotly(plot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 56 rows containing non-finite values (`stat_smooth()`).
#y=mx + b
#y = Sig_th_gsw
#x - decy
#coeffiecients: intercept = b, decy = m
#Sig_th_gsw = -0004*decy + 33.4
#Lab Assignment 1. Pick a question (include hypothesis)
How does oxygen look like at depth and over time? H1 : Oxygen reaches a max in the upper 1000m before decreasing with further depth
o2_h = hydroS_corrected %>%
filter(`O2(1)`!= -999)
o2_h =
ggplot()+
geom_point(aes(x=`O2(1)`,y = Depth))
Potential Questions :
How do temperature, salinity, and sigma-theta co-vary? Is there a relationship betwen sigma-theta and oxygen? Is there a relationship of any of the parameters with depth? With time? Within a depth range over time? Are there seasonal differences in any of the parameters?
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)
## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00
You can also embed plots, for example:
Note that the
echo = FALSE parameter was added to the code
chunk to prevent printing of the R code that generated the plot.